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ABSTRACT 

Many astrophysical flows occur in inhomogeneous (clumpy) media. We 
present results of a numerical study of steady, planar shocks interacting with 
a system of embedded cyhndrical clouds. Our study uses a two-dimensional 
geometry. Our numerical code uses an adaptive mesh refinement allowing us to 
achieve sufficiently high resolution both at the largest and the smallest scales. We 
neglect any radiative losses, heat conduction, and gravitational forces. Detailed 
analysis of the simulations shows that interaction of embedded inhomogeneities 
with the shock/post shock wind depends primarily on the thickness of the cloud 
layer and arrangement of the clouds in the layer. The total cloud mass and 
the total number of individual clouds is not a significant factor. We define two 
classes of cloud distributions: thin and thick layers. We define the critical cloud 
separation along the direction of the flow and perpendicular to it distinguishing 
between the interacting and noninteracting regimes of cloud evolution. Finally 
we discuss mass-loading and mixing in such systems. 
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1. INTRODUCTION 

Mass outflows play a critical role in many astrophysical systems ranging from stars to the 
most distant active galaxies. Virtually all studies of mass outflows to date have focused on 
flows in homogeneous media. However, the typical astrophysical medium is inhomogeneous 
with the "clumps" or "clouds" arising on a variety of scales. These inhomogeneities may 
arise due to initial fluctuations of the ambient mass distribution, the action of instabilities, 
variations in the flow source, etc. Whatever the origin of the clumps their effect can be 
dramatic. The presence of inhomogeneities can introduce not only quantitative but also 
qualitative changes to the overall dynamics of the flow. 

A number of studies have attempted to understand the role of embedded inhomo- 
geneities via (primarily) analytical methods (Hartquist et al. 1986), (Hartquist & Dyson 
1988), (Dyson & Hartquist 1992), (Dyson & Hartquist 1994). In these pioneering works it 
was suggested that interactions of the global flow with inhomogeneities may cause significant 
changes in the physical, dynamical, and even chemical state of the system. Two major con- 
sequences of the presence of clumps are mass-loading (i.e. seeding of material, ablated from 
the surface of inhomogeneities, into the global flow) and transition of the global flow into 
a transonic regime irrespective of the initial conditions. The papers cited above considered 
the potential effects of mass-loading on the global properties of a number objects in which 
inhomogeneities can be resolved. Such objects include planetary nebulae, e.g. NGC 2392 
(O'Dell et al. 1990), (Phillips & Cuesta 1999), and NGC 7293 (Burkert & O'Dell 1998), and 
Wolf-Rayct stars, and primarily RCW58, which is believed to be mass-loading dominated 
(Hartquist et al. 1986). 

A number of numerical studies of single clump interactions have been performed ((Klein 
et al 1994) (hereafter KMC), (Anderson et al. 1994), (Jones et al. 1996), (Gregori et al. 
1999), (Gregori et al. 2000), (Jun & Jones 1999), (Miniati et al. 1999), (Lim & Raga 
1999)). In these papers the basic hydrodynamics or MHD of wind-clump and shock-clump 
physics have been detailed (often with microphysical processes included). A few studies of 
shock waves overrunning over multiple clumps exist in the literature (e.g. (Jun et al. 1996)). 
A detailed study of multiple clumps however, where an attempt is made to articulate basic 
physical processes and differentiate various parameter regimes, has not yet been carried out. 
In this paper, (and those which follow), we address the problem of clumpy flows providing a 
description of the dynamics of multiple dense clouds interacting with a strong, steady, planar 
shock. 

The large parameter space and complexity of the problem require significant compu- 
tational effort. To provide the necessary resolution of the fiow we have used an adaptive 
mesh refinement method. This is a relatively new computational technology and because of 
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this we have chosen to investigate so-called adiabatic flows in which radiative cooling is not 
considered. In this regard our approach is similar to that described by Klein et al. (1994) for 
single clumps and we will utilize their results in understanding our multi-clump simulations. 
We note that preliminary results, appropriate to AGN, were presented in (Poludnenko et al. 
2001). 

The plan of the paper is as follows. In Section 2 we describe the numerical experiments, 
the code used, and the formulation of the problem. In Section 3.1 we consider the general 
properties of the shock-cloud interaction in the context of the multi-cloud systems, primarily 
we focus on the four major phases of the interaction process. In Section 3.2 we discuss the 
role of cloud distribution in determining the dynamics of the system evolution. In Section 
3.3 we define several key parameters, that allow us to distinguish between various regimes of 
shock-cloud interaction. Finally, in Section 3.4 we address the issue of mass-loading in such 
systems. 



2. NUMERICAL EXPERIMENTS 
2.1. Description of the Code Used 

The code we used for this project is the AMRCLAW package which implements an 

adaptive mesh refinement algorithm for the equations of gas dynamics ((Berger & LeVeque 
1997), (Berger & Jameson 1985), (Berger & Colella 1989), (Berger & Oliger 1984)). In the 
AMRCLAW approach, the computational domain is associated with a logically rectangular 
grid that represents the lowest level of refinement (level 1) and that embeds the nested 
sequence of logically rectangular meshes with finer resolution (levels 2,3,...). The temporal 
and spatial steps of all grids at a level L are refined with respect to the level L-1 grids by the 
same factor, typically 4 in our calculations. The mesh ratios Ai/Ax and At/ Ay are then 
the same on all grids, ensuring stability with explicit difference schemes. 

The core of the code - the AMR module - scans each refinement level every k time steps 
and regenerates all nested higher level grids in order to track the moving features of the 
fiow. Two criteria are used to define cells requiring refinement: Richardson extrapolation 
and steepest gradient. The first criterion ensures that the local truncation error does not 
exceed some predefined tolerance. This is done via comparison of the solution obtained by 
taking 2 time steps on the existing grid with one computed by taking 1 time step on a grid 
that is twice as coarse in each direction. The second criterion ensures that the maximum 
of the gradients of all state variables does not exceed some predefined value and guarantees 
that sufficient refinement is achieved in such regions of the flow as shock waves, boundary 
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layers, etc. Flagged cells are then organized into rectangular grid patches in a manner 
that provides a reasonable compromise between the size and the total number of individual 
patches. Finally, the AMR module of the code ensures that global conservation is preserved 
at grid interfaces via introduction of a conservative flux correction. 

After the grid hierarchy is formed, each grid is forwarded to the integration module. This 
module considers every grid as an independent physical domain. The boundary conditions 
are obtained cither from the physical boundary conditions of the computational domain or 
via interpolation from the neighbouring cells of the next lowest refinement level, depending 
on the grid location. This approach allows us to separate logically the AMR and integration 
modules, which facilitates incorporation of new features into the code. The integration 
proceeds by grid level starting at level 1, which is integrated over a time step, then at level 2 
(it should be integrated over Ri—X'^/xi— yi/yi—t'z/ti time steps to catch up), and so on. The 
solution on each grid is advanced via a second-order accurate Godunov-type finite volume 
method in which second-order accuracy is achieved via fiux-limiting and proper consideration 
of transverse wave propagation. The mult i- dimensional wave propagation algorithm is based 
on the traditional dimensional splitting with the Riemann problem solved in each dimension 
by means of a Roe-approximate Riemann solver (LeVeque 1997). It should be noted, that 
our implementation of the Riemann solver, based on the Roe linearization, does not use any 
additional procedures to ensure satisfaction of the entropy condition, as usually employed 
for this type of Riemann solver. Our analysis shows that the numerical diffusion present in 
the system is sufficient to prevent entropy- violating waves from propagating in the system. 

The hydrodynamic equations we solve are appropriate to a single-fluid system, although 
a passive tracer is introduced in order to track advection and mixing of the cloud material. 
This was implemented as an additional wave family in the Roe solver. 

Our numerical experiments were performed on a coarse grid with the resolution of 
50x100 cells and with the maximum number of reflnement levels equal to 3 (meaning that 
the coarse grid associated with the computational domain embeds not more than two nested 
higher resolution levels). Each higher level has a temporal and spatial step refined by the 
factor of 4 in comparison with the next lowest level and we kept this refinement ratio constant 
for all levels. Such setup provides the equivalent resolution^ of 800x1600 cells. In order to 
facilitate comparison of our numerical experiments with those of KMC, we will describe the 
resolution not in terms of the equivalent resolution but in terms of the number of cells that 
fit in the original maximal cloud radius ao, following the convention of KMC. Then all of 



By equivalent resolution hereafter we mean the resolution of a uniform grid covering all of the compu- 
tational domain and possessing the temporal and spatial step of the highest refinement level. 
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the runs described here in our paper have 32 cells per cloud radius. 

KMC suggested that a minimum resolution of 120 cells per cloud radius is necessary. 
We have performed the simulations of the cloud-shock interaction with the resolution of 
120, 75, and 55 cells per cloud radius. Although we will not describe the details of those 
runs in this paper, the principal difference between the cases with maximum and minimum 
resolution, i.e. 120 and 32 cells per cloud radius, is the rate of instability formation at the 
boundary layers^. This does not seem to have any significant effect on the global properties 
of the interaction or the averaged characteristics of the individual cloud ablation processes. 
Therefore, we find the resolution of 30 cells per cloud radius and above to represent accurately 
the global properties of the interaction process under consideration. Moreover, 30 cells per 
cloud radius is a reasonable compromise between maximizing the size of the computational 
domain and capturing as many small-scale features of the interaction process as possible. 

Finally, another aspect of this problem is the connection between the spatial resolution 
(which naturally sets the smallest scale resolvable in the simulations) and the diffusion and 
thermal conduction length scales. As we will see in section 3.1.4, viscous diffusion and 
thermal conduction in a real physical system operate at length scales comparable to the size 
of a computational cell at the highest refinement level used in our simulations. Therefore, in 
a real system, features smaller than the ones that can be resolved with our resolution could 
not survive over the dynamical time scales relevant to the problem. We will address this in 
greater detail when we discuss the mixing phase of cloud evolution. 

2.2. Formulation of the Problem 

We set up a two-dimensional computational volume, associated with the initial condition 

of N different clouds of radius Oj and density pi embedded in the ambient medium of density 
Pa, and an incident shock wave. Since all of the experiments were performed in the Cartesian 
geometry, the clouds are actually cross-sections of the infinitely long cylinders. We will 
address the importance of the cloud shape in more detail in subsequent work where we will 
consider the fully 3-dimensional case of the shock interaction with spherical clouds. Denoting 
the maximum cloud radius present in the system as Qmax, our computational domain is 
25a^a3; X 500^0^;. This allows us to track the dynamical evolution of the system over greater 
temporal and spatial intervals compared to the Qumax x ^^0"max domain, considered by KMC. 



^For the case of lower resolution the lower rate of instability formation may be somewhat compensated 
by the use of the compressive flux limiters. 
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AU our calculations were performed in a fixed reference frame in which both the clouds 
and the ambient medium are stationary at time t = 0. In this reference frame the horizontal 
axis is taken to be the x axis, the vertical axis - y axis. Initially both the clouds and the 
surrounding intercloud medium arc assumed to be in pressure equilibrium and have pressure 
Pq. Typically, the extent of the region, occupied by the cloud distribution at time t = 0, 
is taken to be not more than 30-35% of the horizontal extent of the computational domain 
with Xl offset by 5% from the left boundary of the computational domain and Xji offset 
by 35-40%. Table 1 below, describing the numerical experiments discussed in this paper, 
provides the details of the cloud distribution in each simulation. Figure 1 illustrates the 
setup of the computational domain ai t — 0. 

In the most general case we assume each cloud to have the same nonuniform density 
profile. The clouds have constant density up to a smoothing transition region at the cloud 
edge which is achieved through a linear or tanh(r) function. We typically set the extent 
of the transition region to the outer 20% of a cloud radius and use the tanh(r) - type 
smoothing function. Therefore, the cloud density profile is of the form 

{Pi — Const, < r < Tj 

Pa + Pi , Pa -Pi tanh(r- (1) 
+ ^ tanh(^) ' ^^^^^-^ 

Although there is very little observational data available concerning the internal struc- 
ture of embedded clouds this particular choice of the density profile seems to be a sufficiently 
good approximation to the real physical clouds and inhomogeneities. Burkert and O'Dell 
(Burkert & O'Dell 1998) discussed the evidence for the exponential density profile in the 
cometary knots of NGC 7293 (Hehx nebula) which is similar to the density profile used by 
us. 

In the simple adiabatic interaction of a cloud with a shock wave there are two dimen- 
sionless parameters that completely define the problem: Mach number of the blast wave, 
Ms, and the density contrast between the cloud and the intercloud medium 

X^ ^ (2) 
Pa 

The range of values spanned by the density contrast x can be quite large and is the 
most important parameter of the problem. For the astrophysical situations of interest this 
range can often cover up to 5 orders of magnitude (from 10 to 10^), presenting a significant 
challenge both for the numerical modeling and for the subsequent interpretation and analysis 
of the results. In order to decrease the extent of this dimension of the parameter space, we 
chose a "compromise" value of the parameter x to be 500. Although the runs we discuss 
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in this paper all use this value of the density contrast, we will briefly discuss numerical 
experiments with 10 < x i$ 1000 in the results section, particularly in the context of the 
problem of mass loading. We will provide a more comprehensive study of scaling with density 
contrast in the subsequent work. 

Another important parameter is the shock wave Mach number Mg. We consider a 
planar steady shock wave propagating into the computational domain from the left. Since 
we operate in the reference frame in which both the clouds and the ambient medium are 
stationary, the shock wave Mach number completely defines the shock velocity as well as 
the conditions of the postshock fiow. Using the sound speed of the ambient medium Ca = 
{'yPo/paY^'^, the shock velocity in the stationary reference frame takes the form 

vs-MsC-^Y. (3) 



pa 

Using the Rankine-Hugoniot relations (Landau & Lifshitz 1959) we have the following 
expressions for the postshock conditions^ 



n=(^^^^)»i (7) 



^vs 1 
y + 1 V M 

We assume that the shock wave is strong, so that the condition 

PS — 

PaCl 

is satisfied. 

Shock wave Mach numbers in astrophysical situations can cover a large range of values. 
Fortunately, the problem becomes practically independent on the Mach number for strong 
shocks, i.e. for Ms ~ 10 and above. Indeed, as it can be seen from the shock conditions (4) 
- (6), for Ms < 10 the postshock density pps, pressure Pps, and velocity vps are at most 
within a few percent of their respective values at Ms — > oo. Moreover, recalling that there 
is scale invariance inherent in the hydrodynamic equations under transformations 

t^tMs, v^v/Ms, P^P/Ml (8) 



^In our discussion we assume the perfect gas, i.e. 7 = Const = | for cloud, intercloud, and postshock 
material. 
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the conclusion follows that for Ms ^ 1 the time evolution of a cloud does not depend on the 
Mach number of the shock (Klein et al. 1994)^. Indeed, the results of KMC show that for 
the difference in Ms of 2 orders of magnitude (10 — 1000) time evolution of the system does 
not differ by more than 15%. We will sec that our results fully corroborate the presence of 
Mach scaling in the problem under consideration. 

We assume that the structure of the postshock flow does not change in time for the 
duration of the simulations. An example of such steady postshock flow is the wind from a 
post-AGB star driving a shock with a constant postshock flow structure into a slow wind 
ejected during the previous stages of evolution. This frees us from having to use the pressure 
variation timescale tp to constrain a cloud size, since we can set tp oo. On the other hand, 
for blast waves one cannot assume a steady time independent postshock flow (for example, 
SNR blast waves) and the size of the clouds is constrained by the condition tec <^ tp as 
discussed by KMC (see (Klein et al. 1994)). 

It should be mentioned that the maximum cloud size is still constrained by the condition 
of the shock front planarity. This condition is less restrictive than the one discussed above, 
however it still requires a cloud diameter not to exceed 5 — 10% of the global shock wave 
front radius. This condition is satisfied, for example, in the case of the inhomogeneities, or 
the cometary knots, observed in such planetary nebulae as NGC 2392 and NGC 7293 (e.g. 
(Burkert & O'Dell 1998)). 

The timescale we use to define time intervals in our numerical experiments is the time 
required for the incident shock wave to sweep across an individual cloud, called the shock- 
crossing time, 

tsc = , (yj 

Vs 

where amax = foi' cloud distributions with identical clouds and amax = max{ai) for cloud 
distributions of varying size clouds. 

Due to the scale-invariance of our simulations, one can, using specific values for the 
shock velocity and the size of the inhomogeneities, easily convert the time units used in our 
discussion into the physical ones, tsc is particularly useful to characterize the problem since 
it has clear physical meaning and does not depend on a specific density contrast, which is 
important in the case of systems containing clouds of different density. 

Note that except for the very short period of time when a cloud interacts with the 
shock front, the former finds itself immersed in a post-shock fiow or "wind" the pressure and 



''This conclusion is true with a restriction that the shock speed vs is held fixed. 
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density of which vary only by several percent over the large range of Mach numbers. Since 
KMC showed that the initial interaction with the shock front does not alter the evolution 
of the system for the varying Mach number, the details of the evolution should not change 
after the shock front passed the cloud. Therefore, conclusions about Mach scaling should be 
valid both for the durations of cloud-wind interactions discussed by KMC, and for the much 
longer durations in our experiments. 

One final remark should be made concerning the boundary conditions used in our ex- 
periments. In all runs we imposed a constant inflow at the left boundary, described by the 
postshock conditions, which is determined using the relations (4)-(6), and open boundary 
conditions at the right, top, and bottom boundaries. Those outflow boundary conditions 
were implemented via 0-order extrapolation. 

2.3. Description of the Runs 

All of the runs discussed in this paper contain a Mach 10 shock wave as a part of the 
initial conditions and embedded clouds with the density contrast of 500. Table 1 presents a 
summary of our numerical experiments. 

In addition to the dependence on the shock Mach number and the cloud density contrast 
there are other degrees of freedom present even in the simplest adiabatic case. We considered 
how the dynamical evolution, e.g. rate of momentum transfer from the shock wave and shock 
deceleration, mass loading, mixing of cloud material, etc. of the system depend on 

• the number of clouds present in the system; 

• total cloud mass; 

• spatial arrangement of clouds; 

• individual cloud sizes and masses. 

In most of the runs we constrained ourselves to the case of identical clouds, varying 
only their number and arrangement. Radii of the clouds in all runs except M14^ is 2% of 
the horizontal extent of the computational domain. In order to simplify consideration of the 
dependence on a speciflc cloud arrangement, most runs have a regular cloud distribution, 
where the clouds are placed in the vertices of the mesh, formed by the centers of the clouds in 
the run M14. In addition, we considered a more general case of a random cloud distribution 
with random cloud spatial positions and radii. 
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AU of our numerical experiments were run for about 100 tsc ^- By this time each 
individual cloud has almost completely lost its identity and gained a velocity comparable 
to the velocity of the global flow. Mixing of cloud material with the global ambient flow is 
nearly completed by 100 tsc as well. 

In order to facilitate our analysis, we track temporal evolution of the global averages 
and one-dimensional spatial distributions of several quantities, namely 

• kinetic energy fraction, r^kin = Ekin/Etot, 

• thermal energy fraction, r/jerm = Eterm/Etot, 

• volume filling factor v. 

In order to obtain those quantities from the complex data structure of the adaptive 
mesh simulations, we project the values of the state vector from each grid of the AMR grid 
hierarchy onto a uniform grid with the resolution of the highest refinement level and that 
is associated with the computational domain. Such projection docs not cause loss of data 
or its precision. When this projection is done, we define the global averages of the first two 
quantities above as 



where ^ stands for the quantity being considered, and Ni, Nj are the numbers of cells of the 
projected grid in the x and y direction correspondingly. Such averaging allows us to follow 
momentum transfer from the shock wave to the system of clouds, in the case of {r)kin)2D, and 
heating of the cloud system and intercloud material, in the case of {r]term)2D- 

We also define the one-dimensional spatial averages of those two quantities as 



where ^ again stands for the quantity under consideration. 

Our code follows advection of a passive tracer marking cloud material. In order to 
follow mixing of the cloud material with the global fiow, we define the global average of 



^For comparison, the experiments considered in KMC, that have comparable initial cloud - ambient 
medium density contrast, were run for about 25 tsc- 



N 




(10) 




(11) 
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the volume filling factor {v)^ as the ratio of the total number of cells containing cloud 
material to the total number of cells in the computational domain. We also define the one- 
dimensional spatially averaged volume filling factor {i')id as the variation with the coordinate 
X of the ratio of the number of cells containing cloud material in each vertical row of the 
computational grid to the total number of cells in the vertical dimension. 

3. RESULTS 

3.1. General Properties of the Shock-Cloud Interaction 

Figures 2-5 show the time evolution of a shock wave interacting with a single cloud 
(run Ml), three clouds (run M3), fourteen identical clouds in the regular distribution (run 
MIA), and fourteen clouds of random size in a random distribution (run MlAy). Shown are 
the synthetic Schheren images of the system at four different times for all four sequences. 
Each image is obtained by calculating the density gradient at each point^, plotted on a gray 
scale with the white denoting zero and black - the maximum density gradient. Every image 
in each sequence roughly illustrates transitions between the evolutionary phases discussed 
below. 



3.1.1. Initial Compression Phase 

After initial contact, an external shock transmits an internal forward shock into a cloud. 
This causes cloud compression and heating. At the same time a bow shock forms around the 
cloud. KMC subdivide this phase into two stages: initial transient and shock compression. 
Our numerical experiments show that, in general, their description is applicable for all cloud 
distributions except for the cases when individual clouds are almost in contact at time t = 0. 
However, we do not typically see a reverse shock propagating inside the cloud as they did 
at later stages. The cloud interior seems to be dominated by the forward shock wave which 
prevents a reverse shock from detaching from the back surface of the cloud. The absence 
of the reverse shock is the reason for lower maximum densities in the cloud interior during 
this compression phase compared with KMC: we typically see Pi^max ^ 3pi,o as opposed to 
Pi,max ^ lOpi^Oj quotcd by KMC. Figure 6 illustrates the major flow structures present in the 
system during the initial compression phase. 



^To be more precise, the calculated quantity is the gradient of the density logarithm. This makes the 
images clearer and easier to understand. 
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Propagation of the forward shock in the cloud allows us to define another important 
time scale governing the evolution of the system and defining the duration of the compression 
phase: the cloud crushing time tec This is the time necessary for the internal forward 
shock to cross the cloud and reach its downstream surface 

tec = (12) 
vcs 

In the above expression vcs is the internal forward shock velocity and ttmax is again defined 
as ao in the cases of cloud distributions with identical clouds, and as max{ai) in the cases of 
cloud distributions with clouds of varying size. Following KMC, the velocity of the internal 
forward shock can be written as 

''--^(^)"\fo:F„)"^ (13) 

where vg is the velocity of the external shock. The factor Fgt relates the external postshock 
pressure far upstream with the stagnation pressure at the cloud stagnation point and has 
the form (Klein et al. 1994) 

^ + 1 + 10.7{(7 + I)x}-^/^- ^^^^ 

The factor Fd relates the stagnation pressure with the pressure just behind the internal 
forward shock and has an approximate value of 1.3 determined from numerical experiments 
(Klein et al. 1994). 

While we will primarily use the shock- crossing time as the major time scale, we will 
occasionally give time in terms of the cloud crushing time to facilitate comparison with the 
results discussed by KMC. For this purpose we express the cloud crushing time in terms of 
the shock crossing time. Recalling the definition of tsc (9) we have 

tcc={x"'{^y'\F^iF.,r'^ysc. (15) 
Therefore, for the case of x = 500 

tec = 12 tsc, (16) 
which agrees to about a few percent with the results of the numerical experiments. 



-•^^This was the principal time scale in the study of KMC, although they defined it as as the time necessary 
for the internal forward shock to cross the cloud radius. We have changed the definition in our work since 
the definition of KMC did not actually correspond to the duration of the compression phase. Therefore, tec 
in our work is about twice the tec defined by KMC. 
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The global properties of the flow at this stage are characterized by the onset of individual 
bow shocks around each cloud in a time of order tsc- By the end of the initial compression 
phase those individual bow shocks merge into a single bow shock. 

Finally, the downstream flow, i.e. the flow right behind the external forward shock front, 
is effected by the onset of turbulence in the tails behind the clouds. 



This phase is initiated after the cloud internal forward shock reaches the back of the 
cloud. The two major processes then occur: lateral expansion of the cloud and the onset of 
instabilities at its upstream surface. At this stage Rayleigh-Taylor type instabilities dominate 
at the cloud/ ambient flow interface. These are driven in part by the cloud expansion and 
incipient large-scale fragmentation. The flow downstream with respect to the clouds is 
dominated by Kelvin-Helmholtz instabilities operating in the growing turbulent region. The 
combined action of the lateral expansion and the instabilities causes the clouds to take the 
"umbrella- type" shape and eventually break up. 

In the context of those two processes, the initial cloud separation becomes of key im- 
portance defining the subsequent behaviour of the whole system. We will see below that it 
can be used to distinguish between the two regimes of cloud evolution: interacting and non- 
interacting, and can serve as the basis for classification of cloud distributions. In subsection 
3.3 we will give more rigorous discussion of the role of cloud separation. For now we give a 
qualitative illustration. 



^^It should be noted that a bow wave forms instead of a bow shock if the external postshock flow is 
subsonic, i.e. if 



V Pps / 

With the postshock conditions pps, Pps, and vps determined from the relations (4) - (6), the above criterion 
is satisfied for the following values of the external shock Mach number 



Since in this paper we consider the external shocks, Mach numbers of which are typically above 5.0, we will 
hereafter not consider the possibility of a bow wave formation. 



3.1.2. Re-expansion Phase 





where 



a = 4 - 27(7 - 1), /3 = 7^ - 67 - 7, (5 = 27 -h 2. 
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Clouds, located far enough from each other, are not greatly influenced by their neigh- 
bours and their interaction with the flow proceeds independently as described by KMC. This 
case is illustrated in Figure 7. Compared to the evolution of a single cloud system, shown in 
Figure 2, the two clouds evolve up to the point of their destruction very similarly to the single 
cloud case. However, cloud separations can be small enough for the mutual interaction to 
manifest early during the re-expansion phase, as in Figure 8. This mutual interaction causes 
changes primarily in the flow between clouds. As a result the lateral expansion and growth 
of the Rayleigh- Taylor instabihties in the cloud material is affected. The tails behind the 
clouds are also deformed outwards (see, for example, as well Figure 3). 

The unperturbed supersonic flow that forms behind the external shock wave undergoes 
a transition from a supersonic to a subsonic regime as it passes through a cloud bow shock. 
As a consequence it suffers a signiflcant velocity drop whose magnitude is larger for smaller 
cloud separations due to the larger volume of the stagnation zone in front of the clouds. 
Clouds, acting as the Lavalle nozzles, then cause the flow material to re-accelerate. The 
flow reaches a sonic point next to a cloud core for the regions of the flow adjacent to a 
cloud, and further downstream for the regions of the flow located further from the clouds. 
It is important to note that this re-acceleration results in rarefaction of the flow and a 
gradual decrease both of thermodynamic and dynamical pressure. Eventually, as a result of 
acceleration in the intercloud region, the flow becomes highly supersonic and finally shocks 
down through a stationary shock formed downstream of the clouds to the regime close to 
the unperturbed fiow behind the external shock (see Figures 7-8). 

From the above discussion it is clear that the lateral expansion velocity depends critically 
on the cloud separation. For sufficiently low fiow speeds the cloud material will expand at 
the cloud internal sound speed. With increasing global fiow velocities (or, equivalently, with 
increasing velocities of the external shock front) the lateral expansion velocity will increase 
as well. This velocity is limited, in principle, by the terminal expansion velocity into vacuum. 

For a fixed unperturbed upstream fiow, the fiow velocity near a cloud lateral surface 
(facing the space in between the clouds) will be the highest in the case of a single cloud 
or a cloud located far from the neighbouring ones. With decreasing cloud separation this 
velocity will decrease as well, causing higher dynamical pressure on the lateral surface and, 
therefore, lower lateral expansion velocities. This occurs because the velocity drop across a 
bow shock in the cases of small cloud separations is much larger due to a stronger stagnation 
effect in between the bow shock and the clouds. Therefore, flow adjacent to the cloud does 
not reach velocities as high as in cases of large cloud separations Another way to look at 



^It should be noted that eventually the velocities reached by the flow downstream after passing the region 
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this process is the following. The flow adjacent to the cloud surface passes through a sonic 
point but, in the cases of small cloud separations, densities in the stagnation region are much 
higher. Thus flow densities at the sonic point near the cloud lateral surface are much higher. 
This leads to lower sound speeds and, therefore, lower flow speeds. 

Following KMC, the effective lateral expansion velocity Vexp can be defined as the internal 
cloud sound speed 

„„, = C. = .e,M2fiM, (17) 

7 + i 

where vcs is the velocity of the cloud internal forward shock (13). Our numerical experiments 
prove this to be a very good approximation during almost all of the re-expansion phase. The 
expansion velocity exceeds this value by the end of the re-expansion phase due to stagnation 
pressure in the regions, formed by the Rayleigh- Taylor instability. 

We are now in a position to articulate the temporal evolution of a cloud radius in the 
direction perpendicular to the upstream flow. From the moment of their initial contact 
with the external shock to the moment of their destruction, the clouds first undergo slight 
compression in the direction perpendicular to the flow and subsequently re-expand. KMC's 
analytic model did not explicitly include cloud compression but instead tried to account for its 
effect via a reduced monotonic expansion rate from t — 0. Since a±{t) is intimately related 
to the drag exerted on a cloud by the global flow, the theoretical rate of the momentum 
pickup by a cloud (or the rate of cloud deceleration in the reference frame used by KMC) 
differed from the numerical result. Namely, in Figure 12b of the paper by KMC numerical 
and theoretical results are practically the same up to the time ^ 2.0tcc', when the rate of 
cloud deceleration suddenly increases and the numerical and theoretical results drastically 
diverge. This moment of time corresponds to the beginning of the re-expansion phase, when 
the cloud cross-section starts to increase causing an increase of the rate of the momentum 
transfer from the flow to the cloud. To avoid this problem and simplify an expression for 
a±{t) we use the following form for evolution of a cloud radius normal to the flow, 

\ ao + Cc{t-tcc), tcc<t<tcD. ^ ' 

Here, Cc is given by (17), and tco is the cloud destruction time, defined below in (19). 



between clouds are much higher and, consequently, the strength of the stationary shock downstream is much 
larger in the case of small cloud separations. 
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3.1.3. Cloud Destruction Phase 

Depending on the cloud separation, via the process of re-expansion clouds may come 
into contact and merge into a single coherent structure. This subsequently interacts with 
the flow as a whole and eventually breaks up. Thus for the case of small cloud separations 
we can define the moment of cloud merging as the onset of the cloud destruction phase. For 
large cloud separations in which individual clouds get destroyed before ever merging, it is 
difficult to define the precise onset of the destruction phase as it may be effectively viewed 
as a part of the re-expansion phase. 

We define the end of the cloud destruction phase as the time when the largest cloud 
fragment contains less than 50% of the initial cloud mass. For single cloud systems or systems 
of weakly interacting clouds we define the total time from t — until the end of the cloud 
destruction phase as the cloud destruction time ten, 

tcD = ottcc = ottsc- (19) 

Typically a a; 2.0 consistent with KMC, and using (16) we find a' a; 24 in our simulations. 

In addition to ten there is also a cloud system destruction time tso which we define as 
the time when the largest fragment of a cloud located furthest downstream contains less than 
50% of its initial mass. For thick layer systems (to be described later), including strongly 
interacting cloud distributions, tcD becomes less relevant as a description of the system than 
tsD because tcD < tso- 

3.1.4- Mixing Phase 

After the end of the destruction phase cloud material velocity is still only a small 
fraction of the global flow velocity (see (42) below). The velocity difference promotes Kelvin- 
Helmholtz instabilities at the cloud material - global flow interfaces and, therefore, the 
transition of the system to a turbulent regime. Typically by the beginning of this phase each 
cloud has lost its identity as a result of merging with neighbouring clouds. As the individual 
fragments become smaller and the velocity of the global flow relative to the cloud material 
decreases, Kelvin-Helmholtz instabihties grow faster than the Rayleigh-Taylor type ones. 
This eventually results in complete domination by the former of the small-scale fragmentation 
and causes mixing of cloud material with the fiow (Klein et al. 1994). 

In our numerical experiments, as it can be seen in Figures 2-5, turbulent mixing 
produces a two-phase filamentary system. The appearance of such a two-phase system results 
because our code does not include viscous diffusion or thermal conduction restricting all 



-17- 



dissipative effects to numerical diffusion only. The latter acts at the length scales comparable 
with a cell size at the highest refinement level. 

Real dissipative effects also constrain the overall stability of cold dense plasma embedded 
in a tenuous hotter medium. KMC considered the overall effect of thermal conduction on 
the stabihty of such two-phase media against evaporation. They concluded that the cloud 
ablation time due to evaporation, expressed in terms of the shock-crossing time tsci defined 
in (9) above, has the form 

where F{aQ) is typically of order unity (Klein et al. 1994). Therefore, for the case of our 
simulations, the typical ablation time is tab ~ ^^tsc, or comparable to the cloud destruction 
time. 

One can also estimate an effective depth over which diffusion and thermal conduction 
will disrupt the boundary layer between the two phases over a dynamical time-scale, tsc- 
This can be estimated as follows (see (Kuncic et al. 1996) and references therein). For viscous 
diffusion ^ 

d,.ff - {D,,,ftsct^ - ((^) ^ M^^O ^' (21) 

where Ddiff is the diffusion coefficient and uq is the initial cloud number density. If we 
assume no ~ 1000 cm~^ then, for the cases presented in our simulations, ddiff is about 1% 
of the initial cloud radius, or equivalently is about 1/5 of a cell of the computational domain 
at the highest refinement level. For thermal conduction, the effective depth is 

1 1 i 1 

dterm ~ (D.^rmtsc) ^ = ( ( ^ ) ^ (^) ^ ^ = ^^^^ ^ " ^^2) 

Then dterm is about 6% of ao, or equivalently, about twice the size of a computational cell 
at the highest refinement level. Therefore, should we have included real dissipative effects 
they would destroy the smallest resolvable structures over the dynamically relevant time 
scales. Consequently, any further increase in resolution without providing for the appropriate 
mechanisms, capable to inhibit significantly diffusion and thermal conduction, would not 
provide additional insights into the real physical evolution of a system. 

The importance of the dissipative effects is two-fold. First consider the stability of the 
initial system against destruction due to thermal conduction and diffusion. From the ar- 
guments given above dissipative effects prevent survival of the system for any dynamically 
significant amount of time. As a solution to this problem, KMC suggested that weak mag- 
netic fields inhibit thermal conduction and diffusion. Indeed, as it was shown by (Mac Low 
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et al. 1994), evolution of weakly magnetized clouds during the compression and re-expansion 
phases does not differ significantly from the purely hydrodynamic description. However the 
presence of magnetic fields would raise other issues. During the mixing phase the system 
undergoes transition to turbulence which may amplify the initially dynamically insignificant 
magnetic fields. Turbulence can lower values of the plasma parameter (3 = Pg/Ps ^ 1 or 
even smaller, which may alter the evolution of the system during the later periods of the 
mixing phase. In this respect only a fully magnetohydrodynamic study of the evolution of a 
system of clouds interacting with a strong shock is fully self-consistent (for a series of single 
cloud MHD studies see (Mac Low et al. 1994), (Gregori et al. 1999), (Gregori et al. 2000), 
(Jones et al. 1996), (Miniati et al. 1999), (Lim & Raga 1999), (Jun & Jones 1999)). 



3.2. Role of Cloud Distribution 

In order to characterize the global properties of the shock/cloud system interaction we 
plotted the time evolution of the global quantities defined in section 2.3 for the runs Ml, 
M2, M3, Ab, MIA, and MlAr. Those plots are presented in Figures 9 - 11. 

The important feature of those plots is the striking similarity of the behaviour of systems 
containing similar cloud distributions. The systems containing from one to five clouds ar- 
ranged in a single layer exhibit exactly the same rate of momentum transfer from the global 
flow. This is manifested by the linear rates of fractional kinetic energy {r)kin)2D increase 
from t = up to t = 24:tsc (see Figure 9). The value of the slope for those five cases is 
0.193 ± 1.6%. The thermal energy {rjterm)2D behaves complement arily (see Figure 10). Such 
behavior of single layer systems contrasts that of the multiple layer systems, namely the runs 
Mi4 and Mi4j., which we now discuss. 

The two fourteen cloud runs have different cloud distributions (regular as opposed to 
random), different total cloud mass and different cloud sizes. Nevertheless, the evolution 
of their fractional energies are similar. The rate of the kinetic energy increase during com- 
pression and rc-cxpansion is the same for both M14 and Mi4r and yet is different from that 
single layer cases. The slope in the multi-layer cases is also practically constant throughout 
the two phases with values 0.146 ± 4.5%. 

Note that for all cases the kinetic (thermal) energy reaches its maximum (minimum) 
at the time t = 241 sc, or the time, defined above as the cloud destruction time ten, even 
though for the fourteen cloud runs the cloud system destruction time tsD, defined above in 
subsection 3.1.3, is greater than ten- Thus it seems reasonable to conclude that the cloud 
destruction time is a universal parameter independent on the details of a cloud distribution. 
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After passing through its maximum, the kinetic energy fraction begins to decrease due 
to the transition to turbulence and, consequently, turbulent energy dissipation. It is difficult 
to define a value of the slope for the mixing phase due to the complex nature of the turbulent 
flow but the average rate of kinetic energy dissipation in the system is ^ 0.013 ±27% for the 
single layer systems and ~ 0.015 ± 25% for the multiple layer ones. The proximity of these 
two values (within the standard error) is evidence that the systems have lost any unique 
details of the initial cloud distribution and developed turbulence that depends primarily on 
the rate of energy input at the largest scale, i.e. on the relative velocity of the global flow 
with respect to the cloud material. 

The similarity in behaviour of single vs multiple layer systems is even more prominent 
in the time evolution of volume flUing factors {1^)20 ■ As can be seen in Figure 11, the rate 
of cloud material mixing into the global flow is the same for all single layer systems but is 
different from that in multiple layer ones. The higher mixing rate in the case of multiple layer 
distributions results because upstream clouds pick up momentum faster than the downstream 
ones. Upstream clouds promote destruction of the downstream ones and consequently the 
overall mixing of the system. 

These results lead us to conclude that cloud distribution plays a more important role 
than the number of clouds or the total cloud mass. We use this conclusion in the next section 
as the foundation for classifying possible cloud distributions and defining the general type 
of the cloud system evolution in each category. 

3.3. Critical Density Parameter 

We have seen that the cloud distribution plays the defining role in determining the 
evolution of a shock-cloud system. We now quantify this statement and define criteria for 
determining the behaviour of a given system. 

We define a set of all possible cloud distributions for a given number of clouds N. We 
consider only the clouds of equal or comparable size and density contrast. We define each set 
of cloud distributions Ajv for any given number of clouds to be a set of all possible A^ pairs 
of cloud center coordinates, satisfying two conditions: (1) each pair of clouds is separated 
by some minimum distance Tmin and (2) clouds are confined to a layer extending from the 



^•^The presence of the maximum values in (1^)20 in Figure 11 for all runs is due to the eventual loss of the 
cloud material through the outflow boundaries. 



-20- 



position Xl to the position Xr (see Figure 1): 

ViV > 1 : Ajv = {{xi.yi), l<i <N : = {{xi - x^f + (y^ - y^fY > Tmin ~ 2ao, 

Xi^iXL.XR]]. 

Next, we consider the complete set of all possible cloud distributions A for all possible 
cloud numbers defined as 

A = U {^^}- (24) 

Af>l 

We define within this set A the two subsets: a subset of "thin-layer" cloud distributions A/ 
and a subset of "thick-layer" cloud distributions Am so that 

Ki U Am = A and A/fi Am =0- (25) 

In our numerical experiments those two subsets are associated with the single row and 
multiple row distributions. 

In order to give a precise definition of those two fundamental classes of cloud distribu- 
tions we need to introduce several auxiliary quantities. 



3.3.1. Cloud Velocity and Displacement 

We now estimate the distance that the cloud material will travel before the cloud 
breakup, i.e. within the time ten- 

The equation of motion of a cloud in the stationary reference frame of the unshocked 
ambient medium takes the form 

= ]^CdPps{vps - vcYAc{t), (26) 

where mc is the mass of the cloud, vq is the cloud velocity in the stationary reference frame, 
Cd is the cloud drag coefficient, pps is the undisturbed postshock flow density, and Acit) 
is the cloud cross section area normal to the flow. It should be noted that this equation 
is valid only until the cloud destruction is complete, i.e. until t ~ tco- Prom this point 
on we assume that the drag coefficient Cd ~ 1 which is a rather good approximation for a 
cyhndrical body embedded in a supersonic flow of Mps — 1.31 (see KMC and (Bedogni & 
Di Fazio 1998)). 

Let us assume for a moment that the clouds have flnite extent in the z-direction: zq. 
Note that then 

mc = PoT^alzQ, (27) 
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where po is the cloud density and ao, zq are the cloud dimensions at time t — 0. Moreover, 
the cross section area is 

Ac{t) = 2a^{t)zo, (28) 

where a±{t) is the cloud radius in the direction normal to the flow. Substituting (27) and 
(28) into (26) and using (18) for a±{t) we get the following modified equation of motion 

!^ = ^{«.«-«c}VW. (29) 

at PoT^clq 

This equation describes motion of the cloud as a result of its interaction with the postshock 
wind. However, we also need to account for the velocity that the cloud material acquires 
after its initial contact with the external shock front. This velocity may be comparable to 
the velocity acquired during the compression and re-expansion phases and, therefore, must 
be carefully included into consideration. 

Recall that the initial contact of the incident shock front drives an internal forward 
shock into the cloud with velocity vcs- Cloud material behind the internal shock front gains 
a velocity vc,ps, that can be determined from the Rankine-Hugoniot relations in the usual 
manner. 

Here Mcs is the Mach number of the cloud internal forward shock, which can be expressed 
in terms of the external shock Mach number as follows 

Mcs = ^-Ms( ^V''(F.,F.,)V^ (31) 



Cc,o \ 8/3 

where by Ccfi we denoted the sound speed in the unshocked cloud material and used (13) 
for Vcs- For the simulations discussed in this paper (x = 500 and 7 = 5/3) the internal 
cloud shock Mach number is Mcs — I.86M5 = 18.6. 

Substituting (13) for vcs and (31) for Mcs into (30), and expressing the external shock 
velocity vs in terms of the unperturbed upstream postshock velocity vps by means of (6), 
we obtain the following expression for the velocity of the cloud material due to the cloud 
interaction with the external shock front 



7+1 
8/3 



1/2 



FciFs*)^/^(M|-M(FaF.,)-^) 



vc,PS = VPS ( ^i/2(M|-l) ) = (^2) 



*Note, that from now on we will omit the cloud drag coefficient considering it to be equal to 1. 
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For the case of Mg = 10, F = 0.084. Note that for the hmiting case Ms — > oo the value of 
F remains practically unchanged at 0.083 which corroborates the previously discussed Mach 
scaling (see (8)). 

Finally, making use of the fact that the relative velocity of the postshock flow with 
respect to the cloud is now {vps — vc,ps — vc)-i we can integrate (29) and obtain the following 
form of the cloud velocity 

vc{t) = l ) ^^+(1 r) / (33) 
I " ABit-tccf + Ait-tcc) + c)^ ^ ' ^ 

where we introduced the following quantities 

PoTrao tscX \2\^ + l\ Mp J J ' ^ ' 

1 

R ^ p u/2 37(7-1) y 

C = l2Atsc + (1 - r)-^ (36) 

The unperturbed postshock quantities pps and vps are determined from the conditions (4) 
and (6), Cc is the sound speed in the shocked cloud, defined by (17), and the factor Fgt is 
defined by the relation (14). 

The first quantity A relates the specific momentum of the postshock wind to the cloud 
inertia (mass). Thus it defines the rate of the momentum pickup by a cloud during the 
compression phase, when the cloud dimension transverse to the flow does not increase. The 
second quantity B is the inverse sound crossing time in a compressed cloud, i.e. at the end of 
the compression phase, again for the cloud dimension transverse to the flow. This quantity 
determines the rate of the cloud lateral expansion. Therefore, during the re-expansion phase 
the regular momentum transfer from the wind to the cloud, described by A, is augmented 
by the cloud lateral expansion, described by 5, which comes as an additional factor in the 
quadratic dependence on t. Quantity C ensures continuity of the cloud velocity during the 
transition from the compression to the re-expansion phase. 

Next, integrating (33) from time t = 0.0 up to the cloud destruction time t = ten we 
can determine the displacement of cloud material during the compression and re-expansion 
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phases 

i vps(t-jln(^{l-T)At + l^y 0<t<tcc 

] VPS [t-j (^tan-' ^[ ~tc^f^2c + ^^((^ " ^)^^cc + l))) . tcc<t<tcD 

(37) 

where q = \/ A^BCjA — 1. This allows us to estimate the total displacement a cloud incurs 
before its destruction, called the cloud destruction length, 

LcD = Lc{tcD)- (38) 

In order to get a clearer understanding of the general expressions (33) and (37) let us 
consider the two cases: the case presented in our simulations with Ms = 10 and the limiting 
case of Ms — > oo. We will assume in both cases the density contrast of x = 500.0 and 
7 = 5/3. 

First, we rewrite (33) as 

- f T— «i + «2 ) ) , t <tcc 

Mt) = < > : 2 ^ , . _A (39) 




(^(^^ - 12) bi + + 02) \ tcc<t< tcD 



In the first case of Ms = 10 the coefficients oi, a2, and 61 have the following values 

oi = 1.79 • 10-2; 02 = 1.09; h = 8.35 • 10"^ (40) 

Substituting these into (39) we find that at the end of the compression phase, i.e. at the 
time t = I2tsc, the cloud velocity is 10% of the postshock velocity vps and 7.5% of the 
shock velocity vs- On the other hand, at the end of the re-expansion phase, i.e. at the time 
t = 24:tsc the cloud velocity is 12.66% of vps and 9.4% of vs- 

For the case Ms — > 00 the above coefficients have the values 

ai = 1.83 • 10-2; 02 = 1.09; 61 = 8.51 • 10"^ (41) 

Substitution into (39) gives us the maximum values of the velocity that a cloud can reach in 
the case of an infinitely strong shock: 

vc,max = 10.1 • 10^^fp5 = 7.55 ■ lO-^f^, compressiou pkcise; , . 

Vcmax = 12.8 ■ lO^^fps = 9.57-10^^f5, re — expansion phase. 



^^Note that the assumption here is the same, as in the discussion of Mach scaHng, namely, while increasing 
the shock Mach number, we keep the shock front velocity in the stationary reference frame to be constant. 
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Similarly, we can determine the values of cloud displacement for the two cases, consid- 
ered above. Expression (37) for Lc can be rewritten as follows 

t 1 , / t /ai 



aocA- /n — — +1 , 0<t<tcc 

OqCi C2tan — ■ - C3 , tec <t <tcD 

\tsC V^C4 + C5/ J 

where for the case Ms — 10 the coefficients oi and 02 have the values defined in (40), and 
for Ms 00 the values defined in (41). 

In the case Ms = 10 the coefficients Ci, C2, C3, C4, and C5 have the values 

ci = 1.49; C2 = 104.17; C3 = 10.89; C4 = 9.34 • 10"^; C5 = 113.72. 

Substitution into (43) gives us the displacement that the cloud material undergoes by the 
end of the compression and re-expansion phases: l.Goo and 3.5ao correspondingly. 

In the limiting case Ms — > 00 the values of the coefficients ci, C2, C3, C4, and C5 are the 
following 

ci = 1.5; C2 = 103.22; C3 = 10.9; C4 = 9.43 • 10"^; C5 = 112.56. 

Substituting these coefficients into (43) we find that by the end of the compression phase the 
cloud is displaced by the distance of 1.65ao, whereas by the end of the re-expansion phase 
the displacement is 3.53ao. 

It is clear from the results, obtained above, that both the velocity and cloud displacement 
values in the case Ms — 10 are practically identical to the maximum values, achieved in the 
limiting case of Ms 00. Therefore, our results obtained for the case of a Mach 10 shock 
can be considered as the limiting ones for the cases of strong shocks. 

These results, derived for single clouds or systems with large separation, are in good 
agreement with numerical experiments. Typically the maximum difference between numer- 
ical and analytical values of cloud velocity and position never exceeds 10%. The analytical 
results are usually an overestimate of the numerical ones. This is due to a slight overestima- 
tion of the initial velocity gain after the contact with the external shock front and because 
we assumed the cloud cross-section to be constant during the compression phase, whereas it 
undergoes a small decrease in the experiments. 

Therefore, the maximum distance a cloud can travel before its destruction after the 
initial interaction with a strong shock is 

^CD,max ~ 3.5ao. (44) 
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3.3.2. Critical Cloud Separation 

We first define the average cloud separation, projected on to the direction of the flow 
{Axn) and perpendicular to it (A|/jv), for a given cloud distribution, 



N 

J I N 



(Axjv) - — = (N-m ^ ~ ^^1' ^^^^ 

2 ^ i,j>i 



N 



E 



We can also define a maximum cloud separation projected on to the direction of the flow, or 
the "cloud layer thickness" , 

{AxN)max= HiaX {\Xi-Xj\}. (47) 
i,je{l,N} 



Now we are in a position to give a precise definition of the "thin-layer" and "thick-layer" 
systems. We define a distribution of clouds to belong to the subset A/ if its maximum cloud 
separation {Axi\j)„iax docs not exceed the cloud destruction length Lcd- The distribution 
belongs to the subset in all other cases: 

A/ = {Xn '■ {AxN)max < Lcd}, ^^^-^ 

Am = {Ajv : (Axjv) max > Lcd}- 

A more intuitive way to look at this classification is the following. As we have seen, a 
cloud interacting with the postshock flow re-expands and breaks up before it proceeds into 
the mixing phase. The above criterion tells us if any cloud or a row of clouds will complete its 
destruction phase prior to encountering any other clouds located downstream. The deflnition 
(48) appears to draw rather accurately the line between cloud systems of two types. 

In practice the maximum cloud separation (Axiv)moa; (eq- 47) is simply the thickness 
of the layer of inhomogeneities in a real system and should be compared against the cloud 
destruction length. This thickness can be obtained from the observations of a particular 
object or it can be found analytically, e.g. via consideration of instabilities at the interface 
between two flows. 

Having deflned the two classes, or subsets, of cloud distributions we now consider the 
behaviour of the clouds in each class. First we consider A/, the "single-row" distributions. 
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On average, by the time the clouds are displaced by the distance Lcd-, all of them will be 
destroyed and will proceed to the mixing phase. Thus time of the destruction should be 
approximately too — tsD- 

The question arises whether clouds will interact during the process of re-expansion and 
destruction. We can give a formal criterion for this. Consider two clouds with separation 
(At/Af) = d and {AxN)max ^ Lcd- Both clouds will expand laterally at the velocity Vg^p 
defined in (17). Consequently the time for the clouds to come into contact is 

d — 2(30 / An\ 

Emerge ~ ^ • 
^^exp 

Such re-expansion starts after the cloud compression phase, i.e. after the time tec and cannot 
proceed beyond the cloud destruction time tcD- Therefore, setting tmerge — tcD — tec we 
find the following critical cloud separation transverse to the global fiow 

dcrit ^ 2(ao + Vexp{tcD - tcc))- (50) 

Substituting (17) explicitly for the expansion velocity and (19) for the cloud destruction time 
into (50) we obtain 

tcD-tcc fFciFstY f3^{-f 



In other words clouds whose separation transverse to the fiow is less than dcrit will come 
into contact and merge before their destruction is completed. Therefore, their evolution 
during the destruction phase (and for the most part of the re-expansion phase) can not be 
considered as the evolution of two independent clouds. 

The critical separation does not depend on the global shock Mach number in consistency 
with the Mach scaling, discussed above. Therefore, this parameter is universal for all strong 
shocks and for all possible distributions from the subset A/. For the case 7 = 5/3 and 
X = 500 we find the critical cloud separation to be approximately 

dcrit ~ 4.2ao. (52) 



For cloud distributions from the subset Aj which have an average separation (Ayjv) ^ 
dcrit, the evolution of the system will proceed in the non-interacting regime. On the other 
hand, for the distributions, for which (A|/jv) < da-it, the cloud-cloud interactions are im- 
portant throughout the re-expansion and destruction phases placing them in an interacting 
regime. 



-27- 



It is more difficult to formulate a unified criterion for the behavior of the systems in 
the class A^. When (Axn) > Lcd such systems can be considered as a set of thin layers 
with an average separation Lcd, i-e. each row can be considered as a system from subset 
A/. Consider, for example, the run M14, presented in Figure 4. From Table (1), the average 
separation (Axn) for this run is equal to 7, i.e. (Axn) > Lcd- Indeed, the evolution of the 
leftmost row of clouds proceeds as a simple single row case, and its destruction is completed 
by the time ten- This results in the fractional kinetic energy reaching a maximum at 
time tcD ^ 24t5(7 (see Figure 9). However, it is clear from Figure 4 that the evolution of 
the downstream rows is altered by the destruction of the leftmost one. Therefore, when 
{A.xn) > Lcd one must account for the fact that the destruction of an upstream layer of 
clouds will change the properties of the global fiow for the next, downstream layer. The new 
averaged values of the velocity, density, and pressure in the global flow should then be used 
as an input for the analysis of the downstream cloud layer. 



3.4. Mass loading 

One of the principal questions concerning the effects of shock/cloud-system interactions 
is the role of mass-loading (Hartquist & Dyson 1988). Mass- loading is defined as the feeding 
of material into the global flow by nearly stationary clouds. Analytical studies have predicted 
a number of important changes when mass-loading occurs. The most important of these is 
the transition of the flow to a transonic regime (Hartquist et al. 1986), (Hartquist & Dyson 
1988), (Dyson & Hartquist 1992), (Dyson & Hartquist 1994). In our numerical experiments 
we consider if mass-loading indeed is prominent. 

Mass loading can occur only from time t = up to the moment of cloud destruction at 
time t = tcD- In our experiments the cloud destruction time is fairly short compared with 
the total age of most relevant astrophysical objects. Indeed, cloud destruction is practically 
completed by the time the shock wave reaches the right boundary of the computational 
domain, i.e. by the time the shock wave travels the distance of about 20-30 cloud sizes. 
This could, for example, be compared with clump systems in planetary nebulae. Assuming 
typical size for PNe clouds to be about 100 a.u. (which is the size of cometary knots in 
NGC 7293 (Burkert & O'DeU 1998)), a density contrast 500, and a shock wave velocity 100 
km s~^, we flnd that clouds get completely destroyed within approximately 100 — 150 years. 
This is much less than the typical age of the planetary nebulae (10*^ — 10^ yrs.). 

Thus clouds with low density contrast x-i ~ 10 — 100 can not provide signiflcant mass- 
loading due to the ease in which they are advected and destroyed by the global flow. On the 
other hand, clouds with higher density contrasts Xi > 100 retain their low velocities with 
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respect to the global flow for much longer periods of time and, therefore, may potentially 
be efficient mass-loading sources. However, it should be noted, that this higher relative 
velocity of a cloud increases the efficiency of the instability formation, thereby promoting 
cloud destruction and its mixing with the flow. 

We can also consider the amount of mass seeded into the flow, i.e. stripped of from 
the clouds and assimilated into the global flow, before cloud destruction. Typically, in our 
experiments the amount of seeded cloud material docs not exceed a few percent of the total 
cloud mass, which is unlikely to be enough to switch the flow into a mass-loaded regime. 
Figure 12 shows the distribution of cloud material along the direction of the flow or, to be 
more precise, the distribution of the parameter {i')id{x) for the three cloud run M3. There 
the clouds have the separation (Aun) 0.95dcrit- The first graph corresponds to the end 
of the compression phase, while the second corresponds to the end of the destruction phase. 
The graphs show that cloud material remains localized in the vicinity of the cloud cores 
until the moment of cloud destruction and the system docs not exhibit any significant mass- 
loading. Moreover, the graphs 3 and 4 of Figure 12, showing cloud material distribution early 
in the mixing phase, indicate that even after destruction cloud material remains localized 
within the region of about 8 cloud radii and retains almost the same average velocity with 
respect to the global flow. Only further on in the mixing phase does cloud material spread 
significantly. 

Concluding, we may say that for the cloud density contrast values in the range Xi ~ 
10 — 1000 and practically all values of the global shock wave Mach number, the flows are not 
likely to be subject to mass loading. These flows will be dominated by the mixing of cloud 
material with the global flow that occurs after cloud destruction. Systems with very dense 
clouds Xi >> 1000 may provide sites suitable for mass loading. Future numerical studies 
should be able to confirm this. 



4. CONCLUSIONS 

We have numerically investigated the interaction of a strong, planar shock wave with 
a system of inhomogeneities. These "clumps" are considered to be infinitely long cylinders 
embedded in a tenuous, cold ambient medium. We have assumed constant conditions in 
the global postshock fiow, thereby constraining the maximum size of the clouds only by the 
condition of the shock front planarity. Our results are applicable to strong global shocks 
with Mach numbers 3 < Ms ^ 1000. The range of the applicable cloud/ambient density 
contrast values is 10 — 1000. 
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We considered four major phases of the cloud evolution due to the interaction of the 
global shock and postshock flow with a system of clouds. These are: initial compression 
phase, re-expansion phase, destruction phase, and mixing phase. We describe a simple 
model for the cloud acceleration during the first three phases, i.e. prior to its destruction, 
and derive expressions for the cloud velocity and displacement. The results of that model are 
in excellent agreement with the numerical experiments. The difference in the values of cloud 
velocity and displacement between analytical and numerical results is < 10%. The maximum 
cloud displacement due to its interaction with a strong shock (prior to its destruction) does 
not exceed 3.5 initial maximum cloud radii. The maximum cloud velocity is not more than 
10% of the global shock velocity. 

The principal conclusion of the present work is that the set A of all possible cloud 
distributions can be subdivided into two large subsets A/ and Am- The first subset is 
A/, thin-layer systems. This subset is defined by the condition that the maximum cloud 
separation along the direction of fiow, or the cloud layer thickness, is not greater than the 
cloud destruction length {AxN)max < Lcd- The thick-layer systems Am, are defined by 
the condition {AxN)max > Lcd- The evolution of cloud distributions within each subset 
exhibit striking similarity in behaviour. We conclude that the evolution of a system of 
clouds interacting with a strong shock depends primarily on the total thickness of the cloud 
layer and the cloud distribution in it, as opposed to the total number of clouds or the total 
cloud mass present in the system. The key parameters determining the type of the cloud 
system evolution are therefore the critical cloud separation transverse to the fiow dcru (this 
is also the critical linear cloud density in the layer), and the cloud destruction length Lcd- 

For a given astrophysical situation our results indicate that one might determine, either 
from observations or from theoretical analysis, the thickness of the cloud layer (Ax7v)max- 
This will then determine the class of the given cloud distribution, Aj or Am- For cloud 
distributions from the set A/ with average cloud separation (Ay^) > dcru evolution of 
the clouds during the compression, re-expansion, and destruction phases will proceed in the 
noninteracting regime and the formalism for a single cloud interaction with a shock wave (e.g. 
KMC, (Jones et al. 1996), (Mac Low et al. 1994), (Lim & Raga 1999)) can be used to describe 
the system. On the other hand, if the cloud separation is less than the critical distance, the 
clouds in the layer will merge into a single structure before their destruction is completed. 
Though throughout the compression phase they can still be considered independently of each 
other, their evolution during the re-expansion and destruction phases clearly proceeds in the 
interacting regime. 

When the distribution belongs to the subset Am it is necessary to determine the average 
cloud separation projected onto the direction of the fiow (Ax at), defined by (45) above, and 



-30- 



compare it against Lcd'- if (Axjv) > Lcd evolution of the cloud system can be roughly 
approximated as the evolution of a set of distributions from the subset A/ and the above 
"thin-layer case" analysis applies. If, on the other hand, (Ax at) < Lcd (especially if (Ayjv) < 
dcrit) the system evolution is dominated by cloud interactions and a thin layer formalism is 
inappropriate. 

Finally we have considered the role of mass-loading. Here our principal conclusion is 
that the mass-loading is not significant in the cases of strong shocks interacting with a system 
of inhomogeneities for density contrasts in the range 10 — 1000. In part this is due to short 
survival times of clouds under such conditions, and in part due to the very low mass loss 
rates of the clouds even during the times prior to their destruction. Mass loading may well 
be important in higher density clouds (Dyson & Hartquist 1994). 

The major limitation of our current work is the purely hydrodynamic nature of our 
analysis that does not include any consideration of magnetic fields. As it was discussed in 
section 3.1.4, cold dense inhomogeneities (clouds) embedded in tenuous hotter medium are 
inherently unstable against the dissipative action of diffusion and thermal conduction. This 
evaporates the clouds on the timescales comparable to, or shorter than, the timescales of 
the dynamical evolution of the system. It was suggested that the magnetic fields may play 
a stabilizing role against the action of the dissipative mechanisms. Although weak magnetic 
fields, that are dynamically insignificant up to the moment of cloud destruction, can inhibit 
thermal conduction and diffusion, those magnetic fields may become dynamically important 
due to turbulent amplification during the mixing phase. A fully magnetohydrodynamic 
description of the interaction of a strong shock with a system of clouds will need to be 
carried forward in future works. 

This work was supported in part by the NSF grant AST-9702484 and the Laboratory 
for Laser Energetics under DOE sponsorship. 

The most recent results and animations of the numerical experiments, described above 
and not mentioned in the current paper, can be found at www.pas.rochester.edu/~wma. 
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Table 1. Summary of the Runs Discussed 



Run # of clouds ^ Distribution # of rows x-spacing ^ y-spacing ^ 



Ml 1 regular 1 - - 

M2 2 regular 1 - 4 

A2 2 regular 1 - 12 

M3 3 regular 1-4 

A5 5 regular 1-4 

M14 14 regular 3 7 4 

MUr 14 random 3 3.5^ 3.5^ 



^ Total number of clouds present in the system. 

^ Spacing between the centers of clouds in two different rows, projected onto 
the X-axis, in the units of the maximum cloud radius a^ax 

^ Spacing between the centers of clouds in the same row, projected onto the 
y-axis, in the units of the maximum cloud radius ttmax 

^ Maximum absolute spacing between the cloud centers in any direction. 
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Fig. 1. — Setup of the computational domain. Shown is the setup for the run M14. "x- 
spacing" and "y-spacing" are the parameters used in Table 1 for the description of the runs. 
Note: not drawn to scale. 
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Fig. 2. — Run Ml. Time evolution of a system, containing a single cloud and interacting 
with a Ms = 10 shock wave. Shown are the synthetic Schlieren images of the system at 
times 22.47 tsc, 35.23 tsc, 50.54 tsc, 68.40 tsc- 
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Fig. 3. — Run MS. Time evolution of a system, containing three identical clouds and 
interacting with a Ms = 10 shock wave. Shown are the synthetic Schlieren images of the 
system at times 22.47 tsc, 35.23 tsc: 50.54 tsc: 68.40 tsc- 
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Fig. 4. — Run M14. Time evolution of a system, containing fourteen identical clouds in a 
regular distribution and interacting with a Ms = 10 shock wave. Shown are the synthetic 
Schlieren images of the system at times 22.47 tsc, 35.23 tsc, 50.54 tsc, 69.09 tsc- 
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Fig. 5. — Run M14r. Time evolution of a system, containing fourteen clouds in a random 
distribution and interacting with a Ms = 10 shock wave. Shown are the synthetic Schlieren 
images of the system at times 9.71 tsc, 22.47 tsc, 45.43 tsc, 69.09 tsc- 
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Fig. 6. — Flow structure during the initial compression phase. Shown is the Schlieren image 
of the run Ml at time 5.1 tsc- GS - external global forward shock; RB - external reverse 
bow shock; IF - internal forward shock; IR - internal reverse shock; BF - back flow, caused 
by global forward shock convergence on the symmetry axis; VI - primary vortex sheets, 
caused by regular refieetion of the bow shock; Ml - primary Mach reflected shocks, caused 
by Mach reflection of the global forward shock at the symmetry axis; -S'l - primary Mach 
stem (more precisely, two primary Mach stems); Tl - primary triple points; V2 - secondary 
vortex sheets, caused by the primary Mach reflection of the global forward shock (note the 
two stem bulges formed at the base of the secondary vortex sheets near the symmetry axis); 
M2 - secondary Mach reflected shocks; 5*2 - secondary Mach stems; T2 - secondary triple 
points. 
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Fig. 7. — Run A2. Illustration of the non-interacting regime of cloud evolution: interaction 
of a Ms = 10 shock wave with a system of two identical clouds with the cloud center 
separation of 12.0 ~ 2.86 dcrit- Shown are the synthetic Schlieren images of the system at 
times 22.47 tsc: 35.23 tsc: 50.54 tsc: 68.40 tsc- 
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Fig. 8. — Run M2. Illustration of the interacting regime of cloud evolution: interaction of a 
Ms = 10 shock wave with a system of two identical clouds with the cloud center separation 
of 4.0 ao ~ 0.95 dcrit- Shown are the synthetic Schlieren images of the system at times 
22.47 tsc, 35.23 tgc, 50.54 tsc, 68.40 tsc- 
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Fig. 10. — Time evolution of the global average of the thermal energy fraction {r]term)2D for 
the runs Ml, M2, MS, A5, MU, M14^. 
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Fig. 11. — Time evolution of the global average of the volume fiUing factor {v)^ for the 
runs Ml, M2, M3, Ah, MU, MUr. 
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Fig. 12. — Distribution of cloud material along the horizontal dimension of the computational 
domain for the run M3. Shown are the one-dimensional spatial averages of the volume filling 
factor {p)iD at times 12.26 tsc, 25.02 tsc, 37.78 tsc, 50.54 tsc- 



